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CONVERSION  TABLE 


Conversion  factors  for  U.S.  Customary  to  metric  (SI)  units  of  measurement 


MULTIPLY - 

angstrom 

I.OCO  000  X  E  -10 

Meters  (m) 

atmosphere  (normal) 

1.013  25  X  E»2 

Kilo  pascal  (kPa) 

bar 

1.000  000  XE -2 

Kilo  pascal  (kPa) 

barn 

1.000  000  X  E  -23 

meter2  (m2) 

British  thermal  unit  (thermochemicai) 

1.054  350  XE -3 

joule  (J) 

cal  (thecmochemicai)/cm2 

4.184  000  X  E  -2 

mega  jcuie/m2  (MJ/m2) 

calorie  (thermochemicaj) 

4.184  000 

joule  (J) 

calorie  (thermochemicai/g) 

4.184  000  XE -3 

joule  per  kilogram  (J/kg) 

curies 

3.7C0  000  X  E  -1 

giga  becquerel  (Gbc)* 

aegree  Celsius 

tK-t’K  +  273.15 

degree  keivin  (K) 

degree  (angle) 

1 .745  329  X  E  -2 

radian  (rad) 

aegree  Fahrenheit 

tK«taK  + 459.57/1 .3 

degree  keivin  (K) 

electron  voit 

1.302  19  XE -19 

joule  (J) 

erg 

1.000  000  XE -7 

joule  (J) 

erg/  second 

1.000  000  X  E  -7 

watt  (W) 

foot 

3.048  000  X  E  -1 

meter  (m) 

foot-pound-force 

1.355  818 

joule  (J) 

gallon  (U.S.  liquid) 

3.735  412  X  E  -3 

meter3  (m3) 

inch 

2.540  000  X  E  -2 

meter  (m) 

jerk 

1.000  000  XE -9 

joule  (J) 

joule/ kilogram  (J/kg)  (radiation  dose 
absorbed) 

1.000  000 

gray  (Gy) 

kilotons 

4.183 

terajoules 

kip  (1000  Ibf) 

4.448  222  X  E  *3 

newton  (N) 

kip/inch2  (ksi) 

6.894  757  X  E  -3 

kilo  pascal  (kPa) 

ktap 

1.000  000  XE  *2 

newton-secona/m2  (N-s/m2) 

micron 

1.000  000  XE  -5 

meter  (m) 

mil 

2.540  OOOX  E-5 

meter  (m) 

mile  (international) 

1.609  344  X  E  +3 

meter  (m) 

ounce 

2.334  952  X  E  -2 

kilogram  (kg) 

pound-force  (Ibf  avoirdupois) 

4.448  222 

newton  (N) 

pouna-forca  inch 

1.129  348  X  E  -1 

newton-meter  (N.  m) 

pound-force/inch 

1.751  268  X  E  *2 

newton/meter  (N/m) 

pound-forca/foot2 

4.788  026  X  E  -2 

kilo  pascal  (kPa) 

pound-force/inch2  (psi) 

6.394  757 

kilo  pascal  (kPai 

pouna-mass  (Ibm  avoirdupois) 

4.535  924  X  E  -1 

kilogram  (kg) 

pouna-mass-foot2  (moment  of  inertia) 

4.214  011  X  E  -2 

kilogram-meter  2  (kg.  m2) 

pouna/mass/foot3 

1.061  846  X  E  fl 

kilogram-meter3  (kg/m3) 

rad  (radiation  dose  absorbed) 

1.000  000  XE -2 

gray  (Gy)** 

roentgen 

2.579  760  X  E  -l 

coulomb/kilogram  (C/kg) 

shake 

1.000  000  XE -3 

second  (s) 

slug 

1.459  390  X  E-1 

kilogram  (kg) 

torr  (mm  Hg,  0«C) 

1.333  22  X  E-1 

kilo  pascal  (kPa) 

'  The  becquerel  (8q)  is  the  SI  unit  of  radioaccvity;  1  Bq  »1  event/s. 
"  The  Gray  (Gy)  is  the  SI  unit  of  absorbed  radiation. 


TABLE  OF  CONTENTS 


Section  Page 

PREFACE .  iii 

CONVERSION  TABLE .  iv 

LIST  OF  ILLUSTRATIONS .  vii 

1  INTRODUCTION .  1 

2  NUMERICAL  EXPERIMENTS .  3 

2 . 1  DISPERSION  RELATION .  3 

2 . 2  GRAVITATIONAL  EXCHANGE  RAYLEIGH-TAYLOR 

EXPERIMENTS .  4 

2 . 3  RICHTMYER-MESHKOV  EXPERIMENTS .  1  3 

2.4  DISCUSSION .  17 

3  INTERMEDIATE  ALTITUDE  FIREBALL  INSTABILITIES .  2  1 

4  FIREBALL  TURBULENCE  MODELING .  3  0 

5  SUMMARY  AND  CONCLUSIONS .  3  6 

6  LIST  OF  REFERENCES .  3  8 

Appendices 

A  RAYLEIGH-TAYLOR  EXPERIMENTS  AT  2-CM 

RESOLUTION,  2-ZONE  INITIAL  AMPLITUDE .  4  1 

B  RAYLEIGH-TAYLOR  EXPERIMENTS  AT  1  -CM 

RESOLUTION,  4-ZONE  INITIAL  AMPLITUDE .  4  9 

C  RAYLEIGH-TAYLOR  EXPERIMENTS  AT  1  -CM 

RESOLUTION,  2-ZONE  INITIAL  AMPLITUDE .  5  7 

D  RAYLEIGH-TAYLOR  EXPERIMENTS  AT  0.5  CM 

RESOLUTION,  2-ZONE  INITIAL  AMPLITUDE .  6  5 

E  RICHTMYER-MESHKOV  EXPERIMENTS  AT  2-CM 

RESOLUTION .  7  3 

F  RICHTMYER-MESHKOV  EXPERIMENTS  AT  1  -CM 

RESOLUTION .  7  9 


v 


TABLE  OF  CONTENTS  (Continued) 


Appendices  Page 

G  RICHTMYER-MESHKOV  EXPERIMENTS  AT  0.5  CM 

RESOLUTION .  8  5 

H  FIRST-ORDER  VERSUS  SECOND-ORDER  COMPARISON 

OF  200  KT  AT  50  KM .  9  1 


v  1 


LIST  OF  ILLUSTRATIONS 


Figure  Page 

1  Calculational  geometry  and  setup .  5 

2  Growth  rates  from  Rayleigh-Taylor  calculations .  8 

3  Dimensionless  growth  rates  from  Rayleigh-Taylor 

calculations .  10 

4  Viscosity  from  Rayleigh-Taylor  calculations .  1  1 

5  Reynolds  number  from  Rayleigh-Taylor  calculations .  1  2 

6  Sensible  amplitude  from  Rayleigh-Taylor  calculations .  14 

7  Effect  of  initial  amplitude  on  dimensionless  growth  rates  - 

Rayleigh-Taylor  calculations .  1  5 

8  Growth  rates  from  Richtmyer-Meshkov  calculations  .  1  8 

9  Fireball  evolution  -  200  KT  at  50  km .  2  2 

10  Fireball  evolution  -  multi-MT  at  45  km .  2  3 

1  1  Fireball  evolution  -  1  MT  at  60  km .  2  4 

1  2  Energy  contours  at  90  seconds  -  200  KT  at  50  km  - 

equation-of-state  versus  non-equilibrium  chemistry .  2  6 

1  3  Density  contours  at  90  seconds  -  200  KT  at  50  km  - 

equation-of-state  versus  non-equilibrium  chemistry .  2  7 

1  4  Density  contours  at  60  seconds,  multi-MT  at  80  km, 

first-order  versus  second-order .  2  8 

1  5  Effective  atmospheric  diffusion  coefficient  as  a  function 

of  altitude .  3  1 

1  6  Density  contours  at  90  seconds,  200  KT  at  50  km, 

inviscid  versus  turbulent,  Vy-lO6 .  3  3 

l  7  Density  contours  at  90  seconds,  multi-MT  at  45  km, 

inviscid  versus  turbulent,  Vx~106 .  3  4 

1  8  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0-.04  •  n .  4  2 

1  9  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0.06  •  it .  4  3 

2  0  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  -  0.08  *  n .  4  4 

2  1  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  -  0.10  •  n .  4  5 

2  2  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0.14  •  it .  46 

vii 


LIST  OF  ILLUSTRATIONS  (Continued) 


Figure  Page 

23  2.0-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0.26  •  ji .  4  7 

2  4  1-cm  resolution  -  4-zone  initial  amplitude 

wave  number  =  0.04  •  ji .  5  0 

2  5  1-cm  resolution-  4-zone  initial  amplitude 

wave  number  =  0.06  •  n .  5  1 

2  6  1-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.08  •  n .  5  2 

27  2-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.10  •  n . . .  5  3 

2  8  1-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.14  •  rt .  5  4 

2  9  1-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.26  •  n .  5  5 

3  0  1.0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.04  •  re .  5  8 

3  1  1 .0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.06  •  it .  5  9 

3  2  1.0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.08  •  n .  6  0 

3  3  1.0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.10  •  ji . . .  6  1 

3  4  1.0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.14  •  ji .  6  2 

3  5  1.0-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.26  •  n .  6  3 

3  6  0.5-cm  resolution, -2-zone  initial  amplitude 

wave  number  =  0.04  •  ji .  6  6 

3  7  0.5-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.06  •  ji .  6  7 

3  8  0.5-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.08  •  ji .  6  8 

3  9  0.5-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.10  •  ji .  6  9 


vui 


LIST  OF  ILLUSTRATIONS  (Continued) 


Figure  Page 

40  0.5-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0.14  •  n . . .  7  0 

4 1  0.5-cm  resolution  -  2-zone  initial  amplitude 

wave  number  =  0.26  •  * .  7  1 

42  2-cm  resolution-  2-zone  initial  amplitude 

wave  number  =  0.08  •  « .  7  4 

43  2-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.16  •  it . .  5  5 

44  2-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.20  •  it .  7  6 

45  2-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.36  •  it  — .  7  7 

46  1-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.08  •  it .  8  0 

47  1-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.16  •  it .  8  1 

4  8  1-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.20  •  n .  8  2 

49  1-cm  resolution,  2-zone  initial  amplitude 

wave  number  =  0.36  •  it .  8  3 

50  0.5-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.08  •  it . .  8  6 

5 1  0.5-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.16  •  n .  8  7 

5  2  0.5-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.20*  it .  8  8 

53  0.5-cm  resolution,  4-zone  initial  amplitude 

wave  number  =  0.36  •  it .  8  9 

54  200  KT  at  50  km.  Pressure  contours  at  10  seconds .  9  2 

5  5  200  KT  at  50  km.  Density  contours  at  10  seconds .  9  3 

5  6  200  KT  at  50  km.  Velocity  magnitude  contours  at  10 

seconds . . .  9  4 

5  7  200  KT  at  50  km.  Pressure  contours  at  30  seconds .  9  5 

5  8  200  KT  at  50  km.  Density  contours  at  30  seconds .  9  6 

5  9  200  KT  at  50  km.  Velocity  magnitude  contours  at  30 

seconds .  97 


IX 


LIST  OF  ILLUSTRATIONS  (Continued) 


Figure  Page 

6  0  200  KT  at  50  km.  Pressure  contours  at  60  seconds .  9  8 

6  1  200  KT  at  50  km.  Density  contourst  at  60  seconds .  9  9 

6  2  200  KT  at  50  km.  Velocity  magnitude  contours  at 

60  seconds .  100 

6  3  200  KT  at  50  km.  Pressure  contours  at  60  seconds .  101 

64  200  KT  at  50  km.  Density  contours  at  60  seconds .  102 

6  5  200  KT  at  50  km.  Velocity  magnitude  contours  at 

60  seconds .  103 

6  6  200  KT  at  50  km.  Density  contours  at  60  seconds .  104 


x 


SECTION  1 
INTRODUCTION 


Hydrocodes  have  been  applied  to  problems  of  nuclear  fireball  evolution  for  many  years. 
As  the  codes  and  the  computers  on  which  they  ran  became  more  robust,  and  their  builders 
and  users  more  experienced,  the  calculations  have  become  quite  successful  in  modeling 
observed  fireball  behavior.  With  realistic  initial  conditions,  equations-of-state,  and  ambi¬ 
ent  atmosphere  models,  modem  codes  are  capable  of  reasonable  predictions  of  fireball 
growth  and  rise  rates,  debris  distributions,  and  shock /fireball  interactions.  Given  this  suc¬ 
cess  and  current  interest  in  the  structured  environment  of  intermediate  altitude  fireballs, 
it  was  natural  to  attempt  calculations  at  even  higher  resolution. 

Early  hydrocodes  were  first  order  and  therefore,  diffusive.  As  a  consequence,  their  solu¬ 
tions  were  usually  "smooth".  In  the  quest  for  ever  higher  resolution,  and  to  control 
numerical  diffusion  prior  to  turbulence  modeling,  modern  codes  are  employing  improved 
numerics.  SHARC  (&-CUBED  Hydrodynamic  Advanced  Research  Code)  has  kept  pace. 
Within  the  past  two  years,  its  first  phase  (a  Lagrangian  advancement  of  the  momentum 
and  energy  equations)  was  modified  for  improved  accuracy  and  its  second  phase  (remap  of 
mass,  momentum  and  energy)  was  replaced  with  a  second  order  algorithm.  During  this 
same  period,  a  turbulence  option  employing  a  K-e  model  was  also  incorporated.  Applica¬ 
tion  of  the  improved  models  to  fireballs  has  yielded  some  unexpected  results,  specifically, 
the  development  of  large  scale  structures. 

Simply  stated,  the  second  order  results  differ  considerably  from  previous  first  order 
results  and  are  quite  sensitive  to  seemingly  slight  changes  in  the  implementation  of  the 
numerics.  One  major  consequence  of  the  less  diffusive  numerics  is  that  pressure  and  den¬ 
sity  gradients  are  maintained  to  the  extent  that  conditions  are  favored  for  the  growth  of 
instabilities.  Our  current  turbulence  model  has  not  significantly  altered  the  evolution  of 
the  instabilities. 

Higher  order  calculations  with  SHARC  indicate  that  intermediate  altitude  fireballs  are 
unstable.  This  report  describes  the  characteristics  of  the  instabilities,  the  conditions  neces¬ 
sary  for  their  formation,  and  their  evolution  in  time.  In  the  process  of  interpreting  these 
results  we  completed  a  series  of  calculations  of  known  instabilities  to  characterize  the 
interaction  between  zoning  and  higher  order  differencing.  Time  and  budgetary  restrictions 
have  prevented  the  expenditure  of  the  level  of  effort  necessary  to  fully  explore  the  results. 
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The  numerical  experiments  which  demonstrate  the  response  of  finite  difference  codes 
to  known  instabilities  are  described  in  Section  2.  In  Section  3  we  describe  the  characteristics 
and  evolution  of  our  computed  fireball  instabilities.  Section  4  discusses  the  S-CUBED  tur¬ 
bulence  model  and  the  results  of  its  application  to  fireballs.  Section  5  provides  a  summary, 
conclusions  and  recommendations  for  further  research. 
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SECTION  2 

NUMERICAL  EXPERIMENTS 


The  instabilities  observed  in  the  SHARC  fireball  calculations  developed  in  relatively 
complex  flows  involving  large  density  differences,  possible  effects  due  to  compressibility, 
an  imposed  length  scale  due  to  atmospheric  gradients,  and  curved  streamlines.  To  aid  our 
understanding  of  the  development  and  growth  of  instabilities  we  completed  a  series  of  cal¬ 
culations  of  unstable  flows  in  planar  geometries.  Although  we  gained  considerable  insight 
through  these  efforts,  we  raised  not  a  few  unresolved  questions. 

The  Rayleigh-Taylor  instability  (Chandrasekhar,  1981;  Sharp,  1984)  occurs  in  regions 
where  the  pressure  and  density  gradients  are  of  opposite  sign.  Large  gradients  favor  its 
formation  and  growth.  The  best  known  and  most  studied  example  is  where  a  lighter  fluid 
is  attempting  to  support  a  fluid  of  greater  density  against  gravity.  The  slightest  perturba¬ 
tion  causes  the  surface  to  rapidly  deform  and  the  fluids  to  exchange  places.  Another 
unstable  example  from  gas  dynamics  occurs  when  a  shock  accelerates  a  light  fluid  into  a 
heavier  fluid.  The  resulting  instability  is  known  as  the  Richtmyer-Meshkov  or  shock 
excited  Rayleigh-Taylor  instability  (Richtmyer,  1960;  Sturtevant,  1987).  Examples  of  calcu¬ 
lations  from  both  instabilities  are  presented  in  this  section.  Calculations  of  idealized 
Kelvin-Helmholtz  unstable  flows  were  also  planned  but  not  undertaken  due  to  budgetary 
constraints. 


2.1  DISPERSION  RELATION. 


Consider  two  uniform,  viscous,  incompressible  fluids  in  a  gravitational  field.  The  fluids 
are  separated  by  a  horizontal  interface  that  is  perturbed  in  amplitude  by  a  disturbance  of 
wave  number  k.  The  density  of  the  upper  fluid  is  pi  and  the  density  of  the  lower  fluid  is 
P2-  Linear,  small  amplitude  theory  shows  that  the  disturbance  amplitude  Y  satisfies  the 
equation 


$UnY. 

dt 


(1) 


where  t  is  time  and  n  is  growth  rate.  Chandrasekhar  (1981)  provides  a  thorough  analysis  of 
this  situation  and  a  rather  complex  dispersion  relation  for  the  growth  rate.  Duff  (1962) 
gives  the  following  approximation: 
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n  =  (Agk  +  v2k4)l/2  -  vk2. 


(2) 


where  n  is  the  growth  rate,  A  the  Atwood  number  =  (pi  -  P2)/ (pi  +  P2),  g  is  gravitational 
acceleration,  k  is  the  wave  number  and  v  is  the  kinematic  viscosity.  Equation  (2)  has  the 
same  asymptotic  limits  as  the  expression  by  Chandrasekhar,  namely 


n2  — »  Agk  (k-»0)  n  — »  (k-*») 


(3) 


and  approximately  the  same  maximum  growth  (=nm)  at  the  most  unstable  wave  number 
(=km)-  Expressions  for  nm  and  km  are  given  below: 


nm  = 


(Ag)2/3 

2vV3  ' 


km  = 


1 

2 


(4) 


2.2  GRAVITATIONAL  EXCHANGE  RAYLEIGH-TAYLOR  EXPERIMENTS. 

The  CLOUD  code  (Libersky,  1983)  was  used  to  simulate  the  gravity  driven  rollup  of  the 
unstable  interface  between  a  heavy  fluid  overlaying  a  lighter  one.  CLOUD  is  a  finite  differ¬ 
ence,  incompressible  (anelastic  approximation),  hydrocode  formulated  with  the  stream- 
function  /  vortidty  transport  equations.  Its  main  advantage  over  SHARC  is  that  it  is 
extremely  fast  (5p  s/ (cell-cycle).  Its  use  for  these  calculations  was  many  times  more  eco¬ 
nomical  than  a  similar  set  completed  with  the  fully  compressible  SHARC  code  because  the 
latter’s  time  step  would  have  been  limited  by  sound  speed.  It  was  felt,  however,  that  the 
CLOUD  response  would  be  similar  to  that  from  SHARC  because  they  use  essentially  the 
same  advection  algorithm.  The  assumption  here  is  that  numerical  diffusion  in  the  advec- 
tion  phase  is  the  major  source  of  departure  from  the  inviscid  growth  rate. 

Two  groups  of  calculations,  differing  in  the  way  the  interface  was  perturbed,  were  com¬ 
pleted  with  the  CLOUD  code.  In  the  first  group,  the  interface  was  initialized  as  a  sinusoid 
of  finite  amplitude.  In  the  second,  the  interface  was  initially  flat  and  the  surface  distur¬ 
bance  was  allowed  to  develop  in  time  by  means  of  a  sinusoidal  velocity  perturbation  that 
was  normal  to  the  interface.  Figures  la  and  lb,  respectively,  illustrate  the  geometries. 
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la.  Infinitesimal  initial  disturbance.  lb.  Finite  initial  disturbance. 


lc.  Richtmyer-Meskov 


Figure  1.  Calculational  geometry  and  setup. 
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RIGID  BOUNDARY 


Each  group  included  a  series  of  calculations  at  several  different  resolutions;  three  for 
the  finite  initial  disturbance  and  four  for  the  infinitesimal  initial  disturbance.  Within 
each  series  of  calculations,  the  wave  number  of  the  applied  perturbation  was  varied  from 
7t/25  to  jc/2.  Both  groups  of  calculations  were  completed  in  a  two-dimensional  Cartesian 
mesh  100  by  100  cm  on  a  side.  All  calculations  were  completed  at  an  Atwood  number  of 
0.053  (density  ratio  of  0.9)  and  the  density  of  the  upper  gas  was  1.225  x  10*3  gm/cc.  Tables  1 
summarizes  the  zoning  and  the  initial  interface  conditions. 


Table  1.  Initialization  summary  -  finite  amplitude  disturbance. 


Zoning 

Zone  size 

Initial  amplitude 

50x50 

2.0  cm 

2  zones 

100x100 

1.0  cm 

4  zones 

100  x 100 

1.0  cm 

2  zones 

200  x  200 

0.5  cm 

1  zone 

Table  1.  Initialization  summary  - 

infinitesimal  disturbance. 

Zoning 

Zone  size 

Initial  velocity 

50x50 

2.0  cm 

0.1  cm/s 

76x76 

1.3  cm 

0.1  cm/s 

100x100 

1.0  cm 

0.1  cm/s 

200  x  200 

0.5  cm 

0.1  cm/ s 

During  each  cycle  of  each  calculation,  the  mesh  was  searched  for  the  top  and  bottom 
of  the  crests  in  the  mixing  region  and  this  information,  along  with  the  time,  was  written 
to  a  file.  Subsequent  to  the  calculation,  a  least-squares  fit  of  exponential  form 
(y  =  a  •  exp(n  •  T) )  was  applied  to  the  data  to  determine  the  growth  rate  n.  Equation  2 
was  then  solved  for  the  effective  viscosity. 

The  calculations  that  were  initialized  with  an  infinitesimal  disturbance  did  not  respond 
until  the  surface  had  been  sufficiently  deformed  by  the  perturbation  velocity.  For  these 
calculations,  the  time  scale  was  shifted  by  the  startup  time  before  fitting  and  in  addition. 
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the  "a"  term  from  the  fit  was  interpreted  as  the  minimum  perturbation  amplitude  capable 
of  being  seen  by  the  code. 


Before  presenting  the  growth  rate  results,  it  is  instructive  to  discuss  qualitatively  the 
effects  of  zoning  on  the  evolution  of  the  interface.  Appendices  A  through  D  show  the  time 
evolution  of  the  interface  growth  for  several  wave  numbers  from  calculations  initialized 
with  finite  amplitude  disturbances.  The  reader  is  referred  to  Table  2  for  the  contents  of  the 
appendices. 

Table  2.  Appendix  contents. 

Appendix  Problem  zoning 


A  50  x  50 

B  100  x 100 

C  100  x  100 

D  200  x  200 


-  2  zone  amplitude 

-  1  zone  amplitude 

-  2  zone  amplitude 

-  4  zone  amplitude 


It  is  clear  from  the  plots  in  the  appendices  that,  for  the  same  wave  number  disturbance, 
problems  zoned  differently  evolve  differently.  This  is  contrary  to  the  lore  from  the  days  of 
first  order  codes,  which  proposed  that  continued  doubling  of  resolution  would  lead  to 
convergence. 


Several  processes  contribute  to  the  observed  results.  The  first  and  most  basic  is  that, 
although  the  calculations  were  nominally  initialized  with  the  same  wave  number  distur¬ 
bance,  finite  zoning  resulted  in  the  introduction  of  components  of  higher  frequency.  Dif¬ 
ferences  evolve  because  variations  in  zoning  result  in  variations  in  the  dissipative  and 
dispersive  characteristics  of  the  advection.  The  resulting  differences  in  the  density  field 
then  feed  back  to  the  driving  terms  and  the  solutions  further  diverge. 

The  center-concave  crests  in  the  100  x  100  calculations  at  small  wave  numbers  and  the 
corresponding  center-convex  crests  in  the  200  x  200  zone  calculations  are  particularly  strik¬ 
ing.  Also  of  note  are  the  high  frequency  features  at  edges  of  the  waves  in  the  most  finely 
zoned  calculations  and  the  aliasing  at  the  higher  wave  numbers. 


Figure  2  compares  the  growth  rates  determined  from  the  calculations  with  a  finite  ini¬ 
tial  disturbance  with  those  computed  from  the  calculations  initialized  with  a  flat  interface. 
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The  results  are  remarkably  similar  although  a  dear  trend  indicates  that,  for  a  given  wave 
number  and  zoning,  the  calculations  initialized  with  a  flat  interface  had  the  highest 
growth  rates.  Figure  3  exhibits  the  same  data  in  non-dimensional  form  and  indicates  that, 
in  non-dimensional  space,  the  results  essentially  collapse  to  a  line.  In  these  plots,  the 
growth  rates  were  made  dimensionless  by  dividing  by  the  invisdd  growth  rate  (Agk)1/2. 
Dimensionless  wave  numbers  were  formed  by  expressing  them  as  dx/X  =  2jtdx/k,  where  X 
is  the  wavelength  of  the  imposed  disturbance.  The  dimensionless  wave  number  is 
inversely  proportional  to  the  number  of  zones  across  a  wavdength  and  is  indicative  of 
how  well  the  driving  disturbance  is  resolved. 

The  version  of  CLOUD  used  in  this  series  of  calculations  was  invisdd.  That  is,  no  vis¬ 
cous  terms  were  induded  in  the  governing  equations.  The  non-ideal  response  of  the  code, 
however,  can  be  interpreted  by  means  of  an  apparent  viscosity  through  Equation  2. 

Figure  4  shows  the  apparent  viscosity  plotted  against  dimensionless  wave  number  for  both 
series  of  calculations  and  indicates  that  less  resolved  zoning  results  in  a  correspondingly 
high  apparent  mesh  viscosity.  Figure  5  shows  the  same  information  plotted  in  terms  of  a 
disturbance  Reynolds  number  =  n/vk2  (Yih,  1988). 

Most  disturbing  (confusing)  about  these  results  is  the  high  apparent  viscosity  at  the 
smaller  wave  numbers.  Examination  of  Equation  3  indicates  that  for  small  wave  numbers, 
the  growth  rates  should  approach  the  invisdd  rates.  Chandrasekhar  (1981)  stresses  this 
point  and  states  "viscosity  plays  no  role  among  the  very  long  wavelengths".  Several  fac¬ 
tors,  acting  alone  or  in  concert,  have  been  identified  that  can  contribute  to  the  observed 
disparity  between  theory  and  our  numerical  results. 

The  first  involves  the  possibility  that  equation  2  isn’t  adequate  for  our  analysis  and  that 
perhaps  a  more  complicated  model  such  as  the  diffusion  analysis  of  Duff  (1962)  should  be 
employed.  For  example.  Duff  shows  that  diffusion  decreases  the  amplitude  of  mean  den¬ 
sity.  This  effectively  reduces  the  horizontal  gradients,  and  therefore,  the  driving  mecha¬ 
nism  for  vortidty  generation  (g  dp/dx).  For  poorly  resolved  (in  amplitude)  disturbances, 
this  leads  to  strong  damping.  The  same  mechanism  is  also  active  for  long  wavelengths 
because  the  code  responds  to  a  surface  of  large  radius-of-curvature  as  if  it  were  flat.  For 
sinusoids,  the  largest  radii -of-curvature  occur  at  the  crests. 
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Figure  3.  Dimensionless  growth  rates  from  Rayleigh-Taylor  calculations. 
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Figure  4.  Viscosity  from  Rayleigh-Taylor  calculations. 
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Figure  5.  Reynolds  number  from  Rayleigh-Taylor  calculations. 


A  second  possibility  involves  the  fact  that,  because  inviscid  growth  rates  are  small  for 
small  wave  numbers,  a  systematic  error  in  the  procedures  for  growth  rate  determination 
could  result  in  a  large  apparent  viscosity.  One  source  of  systematic  error  that  was  consid¬ 
ered  pertained  to  the  magnitude  of  the  amplitudes  used  in  the  analysis.  The  validity  of 
linear  theory  is  based  on  the  assumption  that  YK«1,  although  Jacobs  and  Catton  (1988,2) 
observed  that  good  agreement  between  theory  and  experiment  is  obtained  for  values  of 
YK=1.  Our  analysis  used  values  of  YK<rc  to  obtain  sufficient  data  for  fitting  at  large  wave 
numbers.  Growth  rate  comparisons  with  results  obtained  for  YK<1  showed  little  differ¬ 
ence,  however. 

For  the  calculations  initialized  with  a  flat  interface,  the  two  parameter  least-squares  fit 
yields,  in  addition  to  the  growth  rate,  a  number  that  can  be  interpreted  as  the  minimum 
disturbance  amplitude  sensed  by  the  code.  These  results,  shown  in  Figure  6,  indicate  that 
CLOUD  responds  to  signals  that  are  approximately  0.7  of  cell  size.  The  shallow  slope  indi¬ 
cates  that  there  is  a  slight  trend  towards  a  larger  ratio  as  the  signal  becomes  less  resolved. 

The  plots  in  Appendices  B  and  C  show  the  interface  evolution  for  identically  zoned  cal¬ 
culations  that  were  initialized  with  amplitudes  of  two  and  four  zone  heights.  Growth  rate 
results  from  these  calculations  are  displayed  in  Figure  7  and  indicate  that  the  interface  per¬ 
turbations  that  were  initially  better  resolved  were  deformed  at  a  slower  rate.  Equation  1 
(and  perhaps,  one's  intuition)  says  the  opposite. 

Although  we  have  no  clear  explanation,  we  suspect  the  discrepancy  is  caused  by  an 
incomplete  initialization  of  the  problem.  We  failed  to  include  the  initial  velocity  field 
associated  with  the  finite  amplitude  disturbance.  These  comments  also  offer  an  explana¬ 
tion  as  to  why  the  largest  growth  rates  were  noted  in  calculations  in  which  interface  was 
allowed  to  grow  from  an  undisturbed  state. 

2.3  RICHTMYER-MESHKOV  EXPERIMENTS. 

The  SHARC  code  was  used  to  demonstrate  the  response  of  a  fully  compressible  finite- 
difference  hydrocode  to  a  hydrodynamically  unstable  interface.  Figure  1c  shows  the  com¬ 
putational  setup  within  the  two-dimensional  Cartesian  mesh.  In  all  calculations,  a  Mach 
1.25  shock  was  driven  from  a  light  gas  into  a  heavy  gas  through  an  interface  that  was  per¬ 
turbed  in  amplitude  by  a  disturbance  of  wave  number  k.  The  shock,  traveling  from  left  to 
right,  was  initialized  20  zones  to  the  left  of  the  interface  and  was  continuously  fed  from 
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Figure  7.  Effect  of  initial  amplitude  on  dimensionless  growth  rates  - 
Rayleigh-Taylor  calculations. 
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the  left  boundary.  The  Atwood  ratio,  based  on  the  density  behind  the  shock  (P2)  and  the 
unshocked  density  (1.225  x  10'3  g/cc)  to  the  right  of  the  interface  (pi)  was  0.91.  Analogous 
to  the  CLOUD  calculations,  the  zonal  resolution  and  the  wave  number  of  the  applied  dis¬ 
turbance  were  varied  to  determine  their  effects  on  the  growth  rate  of  the  disturbance. 
Table  3  summarizes  the  zoning  and  initial  disturbance  amplitude. 

Table  3.  Initialization  summary  -  Richtmyer-Meshkov  experiments. 


Zoning 

Zone  size 

Initial  amplitude 

76  x  50 

2.0  cm 

2  zones 

150  x 100 

1.0  cm 

2  zones 

300x200 

0.5  cm 

4  zones 

Sturtevant  (1987)  and  Youngs  (1984)  have  presented  analyses  of  the  Richtmyer- 
Meshkov  instability.  Sturtevant's  work  was  experimentally  based;  Youngs  numerical. 

Both  employed  Richtmyer's  (1960)  results  that  showed  "if  the  initial  compression  of  the 
interface  and  of  the  fluids  is  taken  into  account,  the  ultimate  rate  of  growth  of  the  corruga¬ 
tion  agrees,  to  within  5  to  10  percent,  with  that  given  by  the  incompressible  theory." 
Because  the  shock  induced  acceleration  of  the  interface  is  short  lived,  g  =  Ui5t  and  Equa¬ 
tion  1  becomes  for  the  impulsive  case 

££  =  kUiAY  (5) 

dt  1 

where  Ui  is  the  post  shock  mean  velocity  of  the  interface  and  Y'  is  the  post  shock  ampli¬ 
tude  of  the  initial  disturbance.  "Thus  in  the  impulsive  case  the  growth  is  linear  in  time 
rather  than  exponential,  and  acceleration  in  either  direction  leads  to  unbounded  growth." 
(Sturtevant,  1987).  Considerable  ambiguity  exists  for  choosing  Y'  and  Ui  as  there  exist  two 
time  scales  for  the  interaction  of  the  shock  with  the  interface,  one  in  the  fast  gas  and  one  in 
the  slow.  We  used  the  small  compression  model  described  in  Sturtevant  for  the  post 
shock  amplitude 

Yi=1.m.  (6) 

Y0  Us 
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where  Us  is  the  incident  shock  velocity  and  Y0  is  the  initial  perturbation  amplitude.  For 
Ui  we  used  2.2  x  104  cm/s,  a  value  numerically  and  theoretically  determined  for  a  flat 
interface.  No  correction  for  the  Atwood  number  was  applied  because  Sturtevant  indicates 
that  the  post-shock  Atwood  ratio  is  not  much  different  from  the  pre-shock,  even  for  strong 
shocks. 

Appendices  E  through  G  show  the  time  evolution  of  the  interface  growth  for  several 
wave  numbers  from  calculations  zoned  with  2  cm,  1  cm,  and  0.5  cm  resolution.  Clearly, 
the  more  finely-zoned  calculations  do  a  better  job  in  maintaining  gradients.  Little  differ¬ 
ence  is  evident  in  the  amplitudes  of  the  corrugations  in  the  calculations  with  0.5  and 
1.0  cm  resolution;  both  of  which  were  initialized  with  an  amplitude  of  2.0  cm. 

The  calculations  with  2.0  cm  resolution  were  initialized  with  a  disturbance  amplitude 
of  4.0  cm.  For  these,  the  absolute  amplitudes  remained  greater,  especially  in  the  high  den¬ 
sity  fingers,  but  the  growth  rates  were  smaller. 

Growth  rates  were  determined  for  each  calculation  using  procedures  similar  to  those 
described  for  the  Rayleigh-Taylor  experiments.  A  linear  least  squares  fit,  instead  of  the 
previously  used  exponential  fit,  was  employed  in  accordance  with  Equation  5.  Figure  8 
summarizes  the  growth  rate  data  for  the  three  levels  of  resolution  used  in  the  calculations 
in  dimensional  and  non-dimensional  form.  Considerable  scatter  is  evident  and  the  results 
in  non-dimensional  form  do  not  collapse  to  a  line.  Equation  5  indicates  that  the  calcula¬ 
tions  with  the  larger  initial  perturbation  (2  cm  resolution)  should  have  exhibited  the  high¬ 
est  growth  rates.  Because  they  did  not,  we  conclude  that  the  growth  rate  differences  are 
due  to  resolution. 

2.4  DISCUSSION. 

The  methodology  presented  in  Sections  2.2  and  2.3  provides  a  means  for  quantifying 
the  non-ideal  response  characteristics  of  hydrocodes.  A  discussion  of  growth  rate  effects 
alone,  however,  is  incomplete.  In  fact,  growth  rate  errors  are  symptoms  of  more  basic 
errors.  Some  sources  of  these  errors  are  briefly  discussed  in  this  section.  Rood  (1987)  pro¬ 
vides  a  thorough  review  of  the  problems  associated  with  advection  algorithms  and  the 
attempts  that  have  been  made  to  reduce  these  errors.  Much  of  the  material  in  the  next  few 
paragraphs  is  based  on  that  article.  Material  in  quotes  is  taken  verbatim. 
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Numerical  advection  algorithms  are  subject  to  both  dissipation  and  dispersion  errors. 
Dispersion  errors  cause  the  production  of  small  waves  and  results  "from  different  Fourier 
components  of  the  original  distribution  propagating  at  different  phase  speeds."  Dissipa¬ 
tion  errors  exhibit  many  of  the  properties  of  diffusion  and  this  diffusion,  in  some  contexts, 
can  dominate  the  problem.  "Odd  order  schemes  tend  to  be  diffusive  and  even  order 
dispersive." 

Advection  algorithms  are  also  subject  to  aliasing  errors  which  arise  "when  it  is 
attempted  to  resolve  high  wave  number  features  on  a  grid  that  is  too  coarse  to  resolve  the 
features."  Aliasing  errors  cause  the  reflection  of  higher  frequency  energy  into  lower  fre¬ 
quencies;  errors  which  are  evident  in  some  of  the  plots  in  the  appendices  to  this  report. 

The  advection  algorithms  used  in  the  SHARC  and  CLOUD  codes  are  monotone.  That 
is,  they  are  designed  to  prevent  the  generation  of  new  maxima  or  minima,  and  thus 
reduce  dispersion  errors.  However,  Rood  (1987)  quotes  that  "no  linear  scheme  of  second 
order  or  higher  accuracy  can  be  made  free  from  dispersion  errors."  Because  phase  errors 
are  worse  for  poorly  resolved  wavelengths,  "distributions  that  are  rich  in  high  wave 
number  components,  high-order  accurate  differencing  may  reduce  the  accuracy  of  the 
solution." 

With  the  above  comments  in  mind,  "it  can  be  argued  that  the  short  wavelengths 
should  be  selectively  diffused  to  eliminate  those  modes  which  are  not  accurately  mod¬ 
eled."  This  implies  that  for  modes  which  approach  the  resolution  of  the  mesh,  subscale 
turbulence  modeling  must  be  invoked  to  represent  subscale  process.  Rood  (1987)  rails 
physical  diffusion  "an  essential  mechanism  in  any  transport  model."  A  turbulence  model 
was  not  invoked  in  any  of  the  numerical  experiments  discussed  here. 

Sharp  (1984)  states  his  belief  that  a  Rayleigh-Taylor  unstable  interface  is  subject  to  the 
Kelvin-Helmholtz  instability  and  he  poses  the  question  as  to  whether  this  property  might 
eventually  lead  to  the  late  time  self-similarity  and  independence  from  initial  conditions 
noted  by  Youngs  (1984).  The  finely  resolved  calculations  presented  in  Appendix  D  show 
the  formation  of  unstable  regions  at  the  edge  of  the  spikes  that  might  result  from  a  Kelvin- 
Helmholtz  instability.  Future  work  could  address  this  question  specifically.  Another  phys¬ 
ical  source  for  these  features  could  be  the  barodinic  generation  of  vortidty  due  to  density 
gradients  normal  to  the  gravitational  field.  This  mechanism  has  been  noted  previously  by 
Klaassen  and  Clark  (1985)  and  Grabowski  (1989). 
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Close  inspection  of  the  plots  from  the  high  resolution,  long-wavelength  CLOUD  calcu¬ 
lations  reveals  that  fine  scale  structure  was  evident  at  problem  initialization.  This  indi¬ 
cates  that,  if  the  late  time  structure  developed  from  these  features,  some  sort  of  scale 
dependent  filter  (read  diffusion)  would  be  necessary  to  suppress  their  growth.  However,  it 
is  likely  that  both  processes,  the  physical  mechanisms  mentioned  in  the  preceding  para¬ 
graph,  and  the  errors  due  to  advection  and  initialization,  contributed  to  the  development 
of  the  observed  structure.  This  makes  specification  of  the  diffusion  model  more  difficult, 
especially  if  the  response  was  driven  mainly  by  initialization  errors. 

Features  similar  to  the  center  concave  corrugations  evident  in  the  long-  wavelength 
calculations  at  1  cm  resolution  were  noted  in  the  cumulus  cloud  calculations  of  Klaassen 
and  Clark  (1985).  These  authors  also  noted  an  associated  downdraft  that  was  not  evident 
in  the  current  calculations.  As  mentioned  earlier,  we  suspect  the  initialization  procedures 
in  the  current  calculations  introduced  wave  number  components  that  grew  at  differential 
rates. 


20 


SECTION  3 

INTERMEDIATE  ALTITUDE  FIREBALL  INSTABILITIES 


Calculations  with  the  SHARC  code  indicate  that  the  tops  of  intermediate  altitude  fire¬ 
balls  are  unstable.  These  results  contrast  with  prior  first  order  results  which,  except  for  a 
few  instances,  revealed  a  smooth  evolution  in  time.  This  section  reviews  the  develop¬ 
ment  and  evolution  of  the  instabilities  in  calculations  of  buoyant,  transitional,  and  fully 
ballistic  fireballs.  Where  possible,  knowledge  and  experience  gained  from  the  idealized 
numerical  experiments  in  the  previous  section  is  applied. 

The  plots  in  Appendix  H  demonstrate  the  development  of  instabilities  in  a  calculation 
of  200  KT  at  50  km.  First  order  results  are  provided  for  comparison.  Differences  are  readily 
apparent,  most  noticeably  in  the  development  of  structure  at  the  fireball  edge  and  in  the 
definition  of  the  torus  region.  Also  note  the  enhanced  gradients  at  the  fireball  top  and 
side.  Although  the  differences  shown  in  the  appendix  are  specific  to  one  yield  and  height- 
of-burst  combination,  they  are  typical  of  those  seen  in  other  similar  comparisons  involv¬ 
ing  buoyant  and  transitionally  ballistic  fireballs. 

Figure  9  shows  the  computed  time  development  of  the  fireball  from  the  same  calcula¬ 
tion  of  200  KT  at  50  km.  The  secondary  thermal  which  develops  at  the  fireball  top  subse¬ 
quent  to  50  seconds  has  not  been  noted  in  any  prior  calculation  of  this  yield  and  height-of- 
burst  combination.  Figures  10  and  11  present  comparable  plots  for  calculations  of  4  MT  at 
45  km  and  1  MT  at  60  km,  respectively.  The  1  MT  calculation  employed  non-equilibrium 
chemistry.  The  other  two  were  completed  with  an  equilibrium  equation-of-state. 

One  trend  identified  in  viewing  Figures  9  through  11  is  that  the  tangential  length  scale 
(with  respect  to  the  fireball  diameter)  of  the  structures  at  the  fireball  tops  decreases  as  the 
height-of-burst  is  increased.  In  addition,  the  largest  growth  rates  are  experienced  directly  at 
the  fireball  top  for  the  two  lowest  bursts.  For  the  60  km  scenario,  the  largest  growth  rates 
develop  slightly  off  axis.  The  fireball  top  deformities  evident  in  the  plots  are  similar  to 
those  shown  by  Klaassen  and  dark  (1985)  in  their  cumulus  cloud  calculations.  The  latter, 
however,  were  limited  to  thermodynamic  and  moisture  fields  and  were  later  shown 
(Grabowski,  1989)  to  result  from  advection  errors. 

We  speculate  that  the  instabilities  noted  in  the  SHARC  calculations  are  not  numerical 
artifacts  but  result  as  the  code  attempts  to  respond  to  true  hydrodynamically  unstable 
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Figure  10.  Fireball  evolution  -  multi-MT  at  45  km, 
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Figure  11.  Fireball  evolution  - 1  MT  at  60  km. 
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situations.  Higher  order  differencing  steepens  the  gradients  near  the  fireball  edge  which  in 
turn  favors  the  growth  of  instabilities  at  larger  wave  numbers  than  possible  in  first  order 
calculations.  Wavelengths  long  with  respect  to  the  thickness  of  the  gradient  region  see  it  as 
a  discontinuity  and  respond  with  growth  at  a  large  fraction  of  the  invisdd  rate.  Wave¬ 
lengths  that  are  small  with  respect  to  the  gradient  region  grow  correspondingly  slower  or 
are  stable,  depending  on  their  size.  Further  work  is  needed  to  prove  this  premise  but,  on 
the  surface  it  seems  consistent  with  published  theory  (McCartor,  et  al.,  1973)  and  our  own 
numerical  experiments. 

Our  original  chemistry  implementation  in  the  second  order  code  used  an  algorithm 
that  moved  total  mass  in  second  order  and  constituent  masses  in  first  order.  Because 
chemical  energy  moves  with  the  mass,  this  algorithm  is  diffusive  in  both  mass  and  energy. 
Figures  12  and  13  compare  energy  and  density  contours  at  90  seconds  from  equation-of- 
state  and  non-equilibrium  calculations  of  200  KT  and  50  km.  Results  from  the  non-equi¬ 
librium  calculation  are  noticeably  more  diffuse  and  the  secondary  thermal  at  the  fireball 
top  is  much  less  developed.  We  attribute  the  differences  to  the  diffusion  in  the  advection 
algorithm  employed  in  the  chemistry  calculation.  This  suggests  that,  since  diffusion  has 
such  a  large  effect,  physically  real  diffusive  process  like  radiation  and  turbulence  should  be 
included  in  fireball  calculations. 

Current  calculations  of  fully  ballistic  fireballs  are  exhibiting  instabilities  unlike  those 
seen  in  other  calculations.  The  features  described  in  the  preceding  paragraphs  formed  in 
regions  where  a  light  gas  was  pushing  a  heavier  gas;  the  classic  Rayleigh-Taylor  unstable 
situation.  In  the  ballistic  calculations,  no  such  region  forms.  However,  instabilities  do 
develop  in  the  heaved  region.  At  first,  we  suspected  that  responsibility  lay  with  the  unre¬ 
alistic  chemical  rate  processes  related  to  the  constituent  diffusion  noted  above.  We 
addressed  this  by  developing  a  new  algorithm  that  transported  the  constituents,  in  addi¬ 
tion  to  the  total  mass,  in  a  second  order  manner.  The  new  advection  algorithm  had  little 
effect  on  the  development  of  the  instabilities. 

Figure  14  compares  density  contours  at  60  seconds  from  first  order  and  second  order 
material  transport  calculations  of  multi -MT  at  80  km.  The  new  advection  algorithm  was 
used  in  the  second  order  calculation.  Structures  in  the  higher  order  calculation  are  clearly 
visible  at  the  later  time.  Dr.  Bill  Shih  of  PRI  suggested  at  the  1989  AESOP  meeting  in  Santa 
Barbara  that  instabilities  could  develop  in  the  heaved  region  of  intermediate  altitude  fire¬ 
balls.  Dr.  Shih  compared  the  flow  in  the  heaved  region  to  jet  flows  with  an  inflection 


25 


XI ISN30 


kx  joniinv 


»W  3001 I i tv 


27 


Figure  13.  Density  contours  at  90  seconds  -  200  KT  at  50  km  -  equation  of  state  versus 
non-equilibrium  chemistry. 
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Figure  14.  Density  contours  at  60  seconds,  multi-MT  at  80  km,  first-order  versus  second-order 


point  in  the  velocity  profiles.  Such  flows  are  known  to  be  unstable  (Yih,  1977).  Further 
research  would  be  necessary  to  definitively  make  this  connection.  Another  possible  source 
of  instability  for  these  flows  is  what  Book  (1984)  calls  a  convective  instability.  It  develops 
when  the  derivative  of  the  entropy  in  the  direction  of  motion  is  less  than  zero.  It  is  possi¬ 
ble  that  such  situations  could  arise  in  flows  which  involve  complex  chemical  rate  pro¬ 
cesses.  It  is  pertinent  to  note  that  there  is  no  unambiguous  evidence  for  instabilities  in  bal¬ 
listic  fireballs. 

So  far  in  our  discussions,  we  have  ignored  the  real  geometry  of  the  problem.  It  is 
assured  that  real  fireball  instabilities  are  three  dimensional.  They  are  subject  to  wind 
shears,  asymmetric  disassembly  effects,  and  three  dimensional  random  perturbations. 

Two  dimensional  axisymmetric  calculations  do  not  account  for  variations  in  the 
azimuthal  direction  and  real  dynamical  effects  such  as  vortex  stretching.  Consequently, 
questions  arise  as  to  whether  results  from  two  dimensional  calculations  are  realistic.  A 
non-linear  analysis  of  the  Rayleigh-Taylor  instability  by  Jacobs  and  Catton  (1988)  indicates 
that  "axisymmetric  instabilities  grow  faster  than  other  shapes  in  their  respective  geome¬ 
tries".  This  suggests  that  current  calculations,  which  force  axisymmetry,  might  overesti¬ 
mate  the  size  of  the  structured  region,  and  thus  influence  researchers  to  overstate  its 
importance. 

Further  questions  come  to  mind  concerning  the  ability  of  a  finite-difference  code  to  rep¬ 
resent  response  realistically  at  all  wave  numbers.  Our  numerical  experiments  in  Section  2 
indicate  significant  damping  at  wavelengths  not  resolved  over  many  zones.  At  best,  there¬ 
fore,  the  high  wave-number  cutoff  is  governed  by  the  zoning,  as  in  first  order  calculations. 
Concomitant  effects  on  the  energy  cascade  from  large  to  small  scales  are  uncertain,  but  to 
prevent  aliasing,  a  subscale  turbulence  model  is  considered  necessary  for  high  resolution 
calculations. 
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SECTION  4 

FIREBALL  TURBULENCE  MODELING 


The  SHARC  turbulence  model  is  based  on  the  well  known  K-e  model,  (Launder  and 
Spalding,  1972),  which  employs  transport  equations  for  the  turbulent  kinetic  energy  K  and 
its  dissipation  rate  e.  The  specifics  of  the  version  used  in  SHARC  are  described  in  the 
reports  by  Barthel  (1985)  and  Pierce  (1986). 

Two  formulations  of  the  turbulence  model  have  been  incorporated  and  tested  in  the 
SHARC  code.  They  differ  in  the  calculation  of  the  turbulent  eddy  viscosity  (Vt  =  C^k2/e). 

In  the  first,  which  follows  the  derivation  of  Launder  and  Spalding  (1972),  Cp.  is  a  constant 
with  a  value  of  0.09.  In  the  second  (Rodi,  1976,  and  Lakshminarayana,  1986),  Cp  is  a  vari¬ 
able  that  depends  on  the  ratio  of  local  production  of  turbulence  to  its  dissipation.  The 
second  formulation,  called  an  algebraic  Reynolds  stress  model,  reduces  to  the  first  for  equi¬ 
librium  turbulence;  i.e.  production  equals  dissipation.  Both  models  have  been  validated 
against  incompressible  jet  flows,  and  the  variable  Cp  formulation  has  been  successfully 
employed  in  calculations  of  precursed  airblast.  To  date,  however,  results  from  applying 
these  models  to  fireballs  have  been  inconclusive. 

Both  formulations  of  the  model  have  been  applied  to  calculations  of  buoyant  (200  KT  at 
50  km)  and  transition  (multi-MT  at  45  km)  fireballs.  A  strong  sensitivity  to  initial  values 
of  K  and  e  was  found  in  calculations  which  employed  the  constant  Cp  formulation.  If  the 
time  scale  for  dissipation  (K/e)  was  too  large,  unrealistically  high  turbulent  energies  were 
generated  very  rapidly  and  the  calculation  ceased  with  energy  errors.  For  small  K/e,  no 
turbulence  was  generated.  Three  variable  Cp  calculations,  which  had  initial  turbulence 
energy  to  dissipation  rate  ratios  (Kq/Eo)  of  105/1,  105/900,  and  107/107,  were  completed  for 
200  KT  at  50  km.  These,  for  reasons  discussed  in  subsequent  paragraphs,  failed  to  generate 
sufficient  turbulent  energy  to  influence  the  fireball  evolution  noticeably. 

The  variable  Cp  model  is  coded  with  limiters  that  prevent  K  and  e  from  falling  below 
their  initial  values.  In  calculations  where  Kq/Eq  was  107/107,  limiters  kept  both  K  and  e  at 
their  initial  values  and  Cp  remained  near  the  equilibrium  value  of  0.09.  The  net  result 
was  a  near-constant  turbulent  eddy  viscosity  of  approximately  106.  This  appears  to  be  an 
extremely  high  value  until  comparison  is  made  with  data  such  as  that  presented  in  Fig¬ 
ure  15  (Reaction  Rate  Handbook,  1979),  which  shows  the  range  of  Vt  as  a  function  of 
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Figure  15.  Effective  atmospheric  diffusion  coefficient  as  a  function  of  altitude. 


altitude  in  the  ambient  atmosphere.  The  value  of  106  falls  well  within  the  range  of  data  at 
507  km. 

Calculations  of  200  KT  at  50  km  and  multi-MT  at  45  km  were  run  with  an  approxi¬ 
mately  constant  Vj  of  106.  Figures  16  and  17  compare  the  90-second  density  contours  with 
their  invisdd  counterparts.  Discernible  but  insignificant  effects  are  evident.  Comparison 
with  Figure  3  indicates  that  the  numerical  errors  noted  for  the  advection  of  multi- 
component  flows  has  a  larger  effect  than  the  106  viscosity. 

The  SHARC  turbulence  model  accounts  for  turbulence  generation  due  to  the  interac¬ 
tion  of  the  density  and  pressure  fields  through  the  term 

G  =  Crt-^-  VP  •  Vp  (7) 

where  P  is  pressure,  p  is  density,and  Vj  is  the  turbulent  kinematic  (eddy)  viscosity.  Crt  is 
a  constant  (Barthel,  1985),  that  until  recently  was  assigned  the  value  of  4/3.  Recent  experi¬ 
ence  (Pierce,  1989)  suggests  a  smaller  value  might  be  more  appropriate.  We  often  refer  to 
Equation  7  as  the  Rayleigh-Taylor  production  term,  although  "enthalpic  production"  (Issa, 
1980)  might  be  more  generally  applied. 

The  variable  Vj  is  a  function  of  the  local  turbulent  kinetic  energy  K  and  dissipation  rate 
e  through  the  relation  Vt  =  C^k2/p.  In  the  algebraic  Reynolds  stress  model,  C(i  is  propor¬ 
tional  to  the  local  production-dissipation  ratio;  Cy.  is  zero  for  very  low  and  high  values  of 
the  ratio  and  peaks  at  the  equilibrium  value  of  0.09  for  ratios  near  1.  In  the  calculations 
described  in  the  preceding  paragraphs,  the  Rayleigh-Taylor  production  was  included  in  the 
calculation  of  C^.  Consequently,  turbulence  production  in  regions  where  the  Rayleigh- 
Taylor  term  was  large  was  artificially  limited. 

Pierce  (1989)  has  suggested  that  the  variable  C^be  used  only  in  connection  with  the  cal¬ 
culation  of  shear  stresses  and  that  it  be  constant  in  the  calculation  of  the  Rayleigh-Taylor 
term.  We  have  not  yet  applied  this  approach  to  intermediate  altitude  fireballs  because  of 
time  constraints  and  uncertainties  in  the  proper  value  for  Crt  .  An  effort  involving  cali¬ 
bration  of  Crt  remains  to  be  completed.  The  first  step  is  to  find  a  well  documented  exper¬ 
iment  that  can  be  used  as  a  test  calculation. 
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Figure  16.  Density  contours  at  90  seconds,  200  KT  at  50  km,  inviscid  versus  turbulent,  VT-106. 


DENSITY  DENSITY 

- - - - - * - - - 1  160.  t - - - 


**  300111 IV 


m*  30nmTY 


34 


Figure  17.  Density  contours  at  90  seconds,  multi-MT  at  45  km,  inviscid  versus  turbulent,  VT-106 


An  alternate,  and  perhaps  more  general,  approach  would  be  to  derive,  in  a  manner 
analogous  to  that  used  for  the  shear  terms,  algebraic  stress  equations  for  the  Rayleigh- 
Taylor  terms.  This  approach  has  been  applied  with  considerable  success  by  Freeman  (1987) 
in  a  one-dimensional  model.  Implementation  of  this  approach  in  SHARC  would  require 
considerable  work,  due  to  differences  in  formulation  and  application  between  SHARC  and 
the  one-dimensional  model. 
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SECTION  5 

SUMMARY  AND  CONCLUSIONS 


A  series  of  numerical  experiments  of  the  Rayleigh-Taylor  instability  were  completed 
with  the  CLOUD  code.  The  calculations  were  invisdd  and  involved  driving  disturbances 
that  were  nominally  of  a  single  frequency.  Results  indicate  that  numerical  viscosity  pro¬ 
vides  a  major  influence  on  the  evolution  of  the  instability  and  that  the  numerical  viscos¬ 
ity  is  decidedly  zoning  dependent.  Diffusive  effects  are  clearly  evident  in  the  least  resolved 
calculations.  In  addition,  continued  doubling  of  the  resolution  did  not  lead  to  the  conver¬ 
gence  that  is  usually  experienced  with  first  order  codes.  Zoning  dependent  differences  in 
the  dispersive  and  dissipative  characteristic  of  the  advection  algorithm  contributed  to 
these  results,  which  were  exacerbated  by  the  introduction  of  higher  than  desired  frequency 
components  at  initialization.  The  latter  were  a  direct  result  of  attempting  to  initialize  a 
sinusoid  on  a  rectangular  grid. 

Anomalous  results  were  experienced  at  both  the  low  and  high  wave  number  limits 
imposed  by  the  grid.  Low  wave  number  disturbances  experienced  the  greatest  viscosities. 
Consequently,  the  invisdd  growth  rates  predicted  by  theory  at  "small"  wave  numbers 
could  not  be  attained.  Finite  zoning  provided  an  upper  limit  on  the  wave  numbers  that 
could  be  resolved,  let  alone  transported.  Aliasing  effects  were  evident  in  calculations  of 
"high"  frequency  disturbances. 

A  comparable  set  of  calculations  of  the  Richtmyer-Meshkov  instability  were  completed 
with  the  SHARC  code.  Growth  rate  effects  were  observed  that  were  qualitatively  similar  to 
those  noted  for  the  CLOUD  calculations. 

Taken  in  aggregate,  the  numerical  experiments  and  our  fireball  results  emphasize  the 
important  influence  of  zoning  and  numerics  on  hydrodynamic  calculational  results.  Most 
importantly,  they  suggest  that  brute  force  attempts  at  high  resolution  calculations  which 
involve  a  large  number  of  fine  zones,  without  consideration  of  physically  real  diffusive 
processes  such  as  turbulence  and  radiation,  will  likely  lead  to  erroneous  results.  Real  pro¬ 
cesses  which  take  place  on  spatial  scales  of  the  order  of  the  mesh  size  must  be  either  explic¬ 
itly  included  or  parameterized  in  hydrodynamic  calculations.  In  addition,  more  work  is 
warranted  to  further  define  the  effects  of  the  energy  cascade  truncation  which  results  from 
the  high  wave  number  cutoff  imposed  by  the  grid.  The  latter  comment  also  applies  to 
azimuthal  wave  number  limitations  imposed  by  axisymmetric  calculations. 
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Two  turbulence  models,  which  vary  in  the  calculation  of  the  turbulent  eddy  viscosity, 
have  been  implemented  in  SHARC  and  have  been  applied  to  fireball  calculations.  Results 
were  inconclusive.  Shortcomings  in  the  models  have  been  identified  and  work  is  contin¬ 
uing  in  these  areas. 
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Figure  18.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  =  0-.04 
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Figure  19.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  =  0.06  •%. 
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Figure  20.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  -  0.08  •n. 
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Figure  21.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  -  0.10  •%. 
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Figure  22.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  =  0.14  *tc. 
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Figure  23.  2.0-cm  resolution  -  2-zone  initial  amplitude 
wave  number  =  0.26 
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Figure  24.  1-cm  resolution  -  4-zone  initial  amplitude 
wave  number  =  0.04  •%. 
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Figure  25.  1-cm  resolution  -  4-zone  initial  amplitude 
wave  number  =  0.06  *7t. 
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Figure  26.  1-cm  resolution,  4-zone  initial  amplitude, 
wave  number  =  0.08  *tc. 
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Figure  27.  1-cm  resolution,  4-zone  initial  amplitude, 
wave  number  =  0.10  *71. 
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Figure  28.  1-cm  resolution,  4-zone  initial  amplitude, 
wave  number  =  0.14 
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Figure  29.  1-cm  resolution,  4-zone  initial  amplitude, 
wave  number  =  0.26  *Jt. 
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Figure  30.  1.0-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.04  •n. 
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Figure  31.  1.0-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.06  *xc. 
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Rayleigh-Taylor  Experiments  with  CLOUD  Code 
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Figure  32.  1.0-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.08  *7t. 
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Figure  33.  1.0-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.10  •%. 
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Figure  34.  1.0-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.14  •  7t. 
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Figure  36.  0.5-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.04  »n. 
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Rayleigh-Taylor  Experiments  with  CLOUD  Code 
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Figure  37.  0.5-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.06  •rc. 
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Figure  40.  0.5-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.14  •%. 
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Figure  41.  0.5-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.26  *7t. 
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Figure  42.  2-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.08  •«. 
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Figure  43.  2-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.16  #ir. 
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Figure  44.  2-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.20  •%. 
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Figure  45.  2-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.36  *7t. 
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Figure  46.  1-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.08  *7t. 
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Richtmyer-Meshkov  Experiments  with  the  SHARC  Code 
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Figure  49.  1-cm  resolution,  2-zone  initial  amplitude, 
wave  number  =  0.36  •%. 
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APPENDIX  G 

RICHTMYER-MESKHOV  EXPERIMENTS  AT  0.5-CM  RESOLUTION 
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Richtmyer-Meshkov  Experiments  with  the  SHARC  Code 
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Figure  51.  0.5-cm  resolution,  4-zone  initial  amplitude, 
wave  number  =  0.16 
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Richtmyer-Meshkov  Experiments  with  the  SHARC  Code 


90 


APPENDIX  H 

FIRST-ORDER  VERSUS  SECOND-ORDER  COMPARISON  OF  200  KT  AT  50  KM 


This  section  provides  comparison  plots  of  first  order  and  second  order  SHARC  calcula¬ 
tions  of  200  KT  at  50  km.  The  plots  clearly  illustrate  the  magnitude  of  the  differences 
brought  about  by  a  second  order  advection  algorithm. 

By  30  seconds  in  the  second-order  calculation,  the  density  gradients  at  the  fireball  edge 
are  much  steeper  than  those  in  the  first-order  calculation,  and  the  speed  contours  suggest 
the  presence  of  an  instability  near  the  fireball  top.  By  60  seconds,  the  instability  is  evident 
in  the  pressure  field  and  a  secondary  thermal  has  begun  to  form  at  the  fireball  top.  The 
minimum  density  in  the  second-order  calculation  is  smaller  by  more  than  a  factor  of  two 
by  this  time.  By  90  seconds,  the  secondary  thermal  is  well  developed  and  a  bulge  of  approx¬ 
imately  the  same  wavelength  has  begun  to  form  at  the  side.  At  two  minutes,  the  torus  in 
the  second-order  calculation  is  still  very  well-defined  and  apparently  undisturbed  by  the 
unstable  region  above. 
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Figure  54.  200  KT  at  50  km.  Pressure  contours  at  10  seconds. 
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?igure55.  200  KT  at  50  km.  Density  contours  at  10  seconds. 
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Figure  56.  200  KT  at  50  km.  Velocity  magnitude  contours  at  10  seconds. 
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Figure  57.  200  KT  at  50  km.  Pressure  contours  at  30  seconds. 
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Figure  59.  200  KT  at  50  km.  Velocity  magnitude  contours  at  30  seconds. 
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Figure  60.  200  KT  at  50  km.  Pressure  contours  at  60  seconds. 
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Figure  61.  200  KT  at  50  km.  Density  contours  at  60  seconds. 
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Figure  62.  200  KT  at  50  km.  Velocity  magnitude  contours  at  60  seconds. 
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Figure  63.  200  KT  at  50  km.  Pressure  contours  at  90  seconds. 
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Figure  64.  200  KT  at  50  km.  Density  contours  at  90  seconds. 
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Figure  65.  200  KT  at  50  km.  Velocity  magnitude  contours  at  90  seconds. 
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Figure  66.  200  KT  at  50  km.  Density  contours  at  120  seconds. 
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